{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo Particle Simulation Example\n",
    "In this example, we will construct a simple geometry to perform a \"Monte Carlo\" style particle simulation. For simplicity, a random walk will be used rather than actual physics. This example will also perform some basic plotting routines to help visualize the random walk of the particle.\n",
    "\n",
    "## 1 - Create Geometry\n",
    "We will use the same geometry as in the \"Pincell\" Example problem. See that example for more information on how to create a geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using ConstructiveSolidGeometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "top =   Plane(Coord(0.0, 0.0, 150.0),  unitize(Coord(0.0, 0.0, 1.0)),  \"reflective\")\n",
    "bot =   Plane(Coord(0.0, 0.0, -150.0), unitize(Coord(0.0, 0.0, -1.0)), \"reflective\")\n",
    "left =  Plane(Coord(-.63, 0.0, 0.0),   unitize(Coord(-1.0, 0.0, 0.0)), \"reflective\")\n",
    "right = Plane(Coord(0.63, 0.0, 0.0),   unitize(Coord(1.0, 0.0, 0.0)),  \"reflective\")\n",
    "up =    Plane(Coord(0.0, 0.63, 0.0),   unitize(Coord(0.0, 1.0, 0.0)),  \"reflective\")\n",
    "down =  Plane(Coord(0.0, -0.63, 0.0),  unitize(Coord(0.0, -1.0, 0.0)), \"reflective\")\n",
    "clad_outer = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4750)\n",
    "clad_inner = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4180)\n",
    "fuel =       InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4096)\n",
    "cells = Array{Cell}(0)\n",
    "regions = Array{Region}(0)\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(left, -1))\n",
    "push!(regions, Region(right, -1))\n",
    "push!(regions, Region(up, -1))\n",
    "push!(regions, Region(down, -1))\n",
    "push!(regions, Region(clad_outer, 1))\n",
    "ex = :(1 ^ 2 ^ 3 ^ 4 ^ 5 ^ 6 ^ 7)\n",
    "push!(cells, Cell(regions, ex))\n",
    "regions = Array{Region}(0)\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(clad_outer, -1))\n",
    "push!(regions, Region(clad_inner, 1))\n",
    "ex = :(1 ^ 2 ^ 3 ^ 4)\n",
    "push!(cells, Cell(regions, ex))\n",
    "regions = Array{Region}(0)\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(clad_inner, -1))\n",
    "push!(regions, Region(fuel, 1))\n",
    "ex = :(1 ^ 2 ^ 3 ^ 4)\n",
    "push!(cells, Cell(regions, ex))\n",
    "regions = Array{Region}(0)\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(fuel, -1))\n",
    "ex = :(1 ^ 2 ^ 3)\n",
    "push!(cells, Cell(regions, ex))\n",
    "bounding_box = Box(Coord(-.63, -.63, -150), Coord(.63, .63, 150))\n",
    "geometry = Geometry(cells, bounding_box);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 - Define Random Walk Routine\n",
    "This routine performs a straw man random walk step by sampling a random distance for the particle to travel before changing direction. In a real application, this would of course have a lot of physics associated with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "move_particle (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function move_particle(ray::Ray, geometry::Geometry)\n",
    "    # Sample distance to travel\n",
    "    distance_to_scatter = rand() * 0.75\n",
    "    \n",
    "    # Determine next intersection distance\n",
    "    new_ray, id, boundary_type = find_intersection(ray, geometry)\n",
    "\n",
    "    # Compute distance travelled by the ray\n",
    "    distance_to_boundary = magnitude( new_ray.origin - ray.origin )\n",
    "    \n",
    "    # Check if scatter happened before boundary crossing\n",
    "    if distance_to_scatter < distance_to_boundary\n",
    "        # Sample new ray direction\n",
    "        new_ray = generate_random_ray(geometry.bounding_box)\n",
    "        \n",
    "        # Move ray to correct location\n",
    "        new_ray.origin = ray.origin + ray.direction * distance_to_scatter\n",
    "    end\n",
    "    \n",
    "    return new_ray\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - Plot Random Walk of Single Particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Plot Background Geometry\n",
    "plot_geometry_2D(geometry, Box(Coord(-0.63, -0.63, 0), Coord(.63, 0.63, 0)), 1000)\n",
    "\n",
    "# Sample new particle\n",
    "ray = generate_random_ray(geometry.bounding_box)\n",
    "old_ray = generate_random_ray(geometry.bounding_box)\n",
    "\n",
    "# Plot particle steps\n",
    "for i=1:50\n",
    "    ray = move_particle(old_ray, geometry)\n",
    "    x = [old_ray.origin.x, ray.origin.x]\n",
    "    y = [old_ray.origin.y, ray.origin.y]\n",
    "    old_ray.origin = ray.origin\n",
    "    old_ray.direction = ray.direction\n",
    "    plot!(x, y, c=:black)\n",
    "end\n",
    "\n",
    "plot!()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 - Plot paths of 100 particles, colored by cell id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(leg=:false, aspect_ratio=1)\n",
    "\n",
    "# Particle Loop\n",
    "for i = 1:100\n",
    "    \n",
    "    # Randomly sample a new particle\n",
    "    ray = generate_random_ray(geometry.bounding_box)\n",
    "    old_ray = generate_random_ray(geometry.bounding_box)\n",
    "    \n",
    "    # Transport Loop\n",
    "    for j=1:50\n",
    "        \n",
    "        # Move particle forward\n",
    "        ray = move_particle(old_ray, geometry)\n",
    "        \n",
    "        # Generate lines to plot\n",
    "        x = [old_ray.origin.x, ray.origin.x]\n",
    "        y = [old_ray.origin.y, ray.origin.y]\n",
    "        \n",
    "        # Determine which cell the particle was travelling in\n",
    "        id = find_cell_id(old_ray.origin, geometry)\n",
    "        \n",
    "        # Update particle\n",
    "        old_ray.origin = ray.origin\n",
    "        old_ray.direction = ray.direction\n",
    "        \n",
    "        # Plot line\n",
    "        if id == 1\n",
    "            plot!(x, y, c=:blue)\n",
    "        elseif id == 2\n",
    "            plot!(x, y, c=:grey)\n",
    "        elseif id == 3\n",
    "            plot!(x, y, c=:white)\n",
    "        else\n",
    "            plot!(x, y, c=:green)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "plot!()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
